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SECTION  I:  GENERAL  PLASMA  THEORY  AND  SIMULATION 


A.  Vortex  formation  and  particle  transport  in  a  cross-field  plasma  sheath 

K,  Theilhaber 

A  final  report  on  thei  subject,  covering  about  two  years  work,  has  been  issued  as  ERL  Memo  No. 
UCB/ERL  M88/21,  with  the  above  title.  The  abstract  follows. 


Abstract 

The  time-dependent  behavior  of  a  transversely  magnetized,  two- 
dimensional  plasma-wall  sheath  has  been  studied  through  particle 
simulations,  with  the  aim  of  modelling  plasma  behavior  in  the  vicinity 
of  the  limiters  and  walls  of  magnetized  plasma  devices.  The  simula¬ 
tions  have  shown  that  the  cross-field  sheath  between  a  wall  and  a 
plasma  is  a  turbulent  boundary  layer,  with  strong  potential  fluctua¬ 
tions  and  anomalous  particle  transport.  The  driving  mechanism  for 
this  turbulence  is  the  Kelvin-Helmholtz  instability,  which  arises  from 
the  sheared  particle  drifts  created  near  the  wail.  Provided  it  is  re¬ 
plenished  by  an  internal  flux  of  particles,  the  sheath  maintains  itself 
in  a  dynamic  equilibrium,  in  which  the  linear  edge  instability,  the 
nonlinear  dynamics  of  the  particles  and  the  outward  particle  diffusion 
all  balance  each  other.  The  sheath  assumes  an  equilibrium  thick¬ 
ness  of  order  lx  ~  5  p,-,  and  maintains  large,  long-lived  vortices,  with 
amplitudes  6<p  ~  2 T,/e,  which  drift  parallel  to  the  wall  at  roughly 
half  the  ion  thermal  velocity.  The  sheath  also  maintains  a  large, 
spatially- averaged  potential  drop  from  the  wall  to  the  plasma,  with 
A<p  «  —  1.5T;/e,  in  sharp  distinction  with  the  unmagnetized  sheath, 
where  the  plasma  potential  is  higher  than  at  the  wall.  Accompanying 
the  long- wavelength  vortices  is  a  spectrum  of  shorter- wavelength  fluc¬ 
tuations,  which  extend  to  |k|  p,-  ~  1  and  ~  vCT,  and  which  induce  an 
anomalous  cross-field  transport.  A  central  result  is  that  the  anoma¬ 
lous  transport  scales  like  Bohm  diffusion,  at  least  when  u/p,  >  2 uet.  At 
lower  densities,  Wp,  <  2U&,  the  diffusion  coefficient  has  an  additional 
factor,  proportional  to  the  density.  These  results  enable  us  to  model 
the  cross-field  sheath  by  a  simple  boundary  condition,  which  relates 
the  particle  flux  through  the  sheath  to  the  edge  density  and  which 
can  be  used  as  input  in  any  model  designed  to  obtain  the  bulk  plasma 
properties. 
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B.  On  Magnetospheric  Hydromagnetic  Waves 
Excited  by  Energetic  Ring-Current  Particles 


Liu  Chen 

Department  of  Electrical  Engineering  and  Computer  Sciences, 
University  of  California  at  Berkeley,  Berkeley,  Calif.  94720 

and 

Akira  Hasegawa 
AT&T  Bell  Laboratories 
Murray  Hill,  New  Jersey  07974 


Abstract.  A  kinetic  theory,  allowing  finite  ion  Larmor  radii,  general  magnetic 
field  geometries  and  plasma  equilibria,  has  been  developed  for  hydromagnetic 
Alfven  waves  excited  within  the  Earth’s  magnetosphere  by  the  storm-time 
energetic  ring-current  particles.  In  particular,  assuming  adiabatically  injected 
particles,  it  is  found  that  the  excitation  mechanism  is  due  to  the  magnetic 
drift-bounce  resonances  of  highly  energetic  protons  and  the  various  predicted 
instability  properties  are  consistent  with  satellite  observations. 


*  Permanent  address:  Plasma  Physics  Laboratory,  Princeton  University,  Princeton,  N.J. 
08544 


This  is  work  done  here  by  Dr.  Lin  Chen  of  PPPL  while  he  was  a  visiting  Miller  Research  Professor, 
Fall  1987. 


INTRODUCTION 


Recent  multiple-satellite  observations  [Takahashi  et  al.,  1985]  of  dayside 
compressional  Pc5  pulsations  have  stimulated  renewed  theoretical  interest  in 
the  excitation  mechanisms  of  magnetospheric  hydromagnetic  waves  during 
geomagnetic  storms.  Briefly  stated,  the  observations  have  the  following 
properties:  (1)  occurrences  mostly  during  the  recovery  phases  of  geomagnetic 
storms,  (2)  large  azimuthal  (East-West)  wave  numbers;  typically,  m~0(l0*), 

(3)  westward  phase  velocities  similar  to  the  ion  magnetic  drift  velocities,  and 

(4)  frequencies  roughly  corresponding  to  those  of  standing  shear  Alfve'n  waves. 
Properties  similar  to  (2),  (3),  and  (4)  were  also  observed  in  storm-time  Pc5 
pulsations  by  STARE  radar  [Allan  et  al.,  1982).  In  one  event,  antisymmetric 
structures  along  the  magnetic  fieid  lines  were  also  identified  [Takahashi  and 
Higbie,  1986;  Takahashi  et  al.,  1987], 

Since  3  {3**  plasma  pressure/magnetic  pressure)  values  are  typically  low. 
0(0.3),  during  the  recovery  phases,  these  observations  appear  to  indicate 
that  the  drift-mirror  instability  mechanism,  as  originally  proposed  by  Hasegawa 
[1969],  and  its  variations  [Walker  et  al„  1982:  Pokhotelov  et  al..  1986]  may  be 
inadequate  in  explaining  these  events.  On  the  other  hand,  excitations  via 
wave-particle  resonances  •  remain  viable  even  at  low  3  values.  Previous 
theoretical  investigations  in  this  direction  are,  however,  limited  to  either  local 
WKB  ansatz  (i.e.,  neglecting  eigenmode  field-aligned  structures)  [Meerson  et  al.. 
1978]  or  the  long-wavelength  drift-kinetic  approximation  [Southwood.  1976]. 
Since  m~0(l0')»l,  the  drift-kinetic  approximation,  which  assumes 


|kj.p|  <  1  (kj.  and  p  being,  respectively,  the  wave  number  perpendicular  to  B 
and  the  Larmor  radius)  becomes,  strictly  speaking,  invalid  for  highly  energetic 
protons.  Furthermore,  the  |kj.p|<l  approximation  also  implies  that  the 
magnetic  drift  frequencies  are  much  smaller  than  the  bounce  frequencies  and. 
hence,  precludes  the  potentially  important  wa”e-particle  interactions  via 
magnetic  drift-bounce  resonances  [Southwood,  1976].  In  the  present  work,  we 
reexamine  this  subject  employing  the  more  consistent  gyrokinetic  equations 
[Antonsen  and  Lane,  1980;  Catto,  et  al.,  1981],  which  includes  effects  of  finite 
ion  Larmor  radii,  general  magnetic  field  geometries,  and  plasma  equilibria  as 
well  as  field-aligned  eigenmode  structures.  In  order  to  isolate  the  responsible 
physical  mechanisms,  we  shall  ignore  pressure  anisotropy  (i.e.,  the  drift-mirror 
mode)  and  concentrate  on  instabilities  excited  via  wave-particle  interactions. 

THEORETICAL  FORMULATION  AND  GYROKINETIC  EQUATIONS 

To  include  finite  ion  Larmor  radius  effects  in  the  formulation  of  low- 
frequency  hydromagnetic  waves,  we  adopt  the  following  gyrokinetic  orderings 
w/fl'-  kuP ■ —  ku/kx—  p/L0 « 1.  and  k±p~0(l);  where  fi  and  p  are. 
respectively,  cyclotron  frequency  and  radius,  and  L0  is  the  equilibrium  scale 
length.  Thus,  within  these  orderings,  the  parallel  wave  length  is  0(L0),  and  one 
must  retain  the  parallel  eigenmode  structures.  Meanwhile,  the  linearized  Vlasov 
equation  can  be  reduced  by  first  transforming  to  the  guiding-center  phase  space 
[X=  x  +  v  X  en/n,  V  =  (e  =  v2/2,  m  =*  VI/2B,  a=»  gyrophase  angle)]  and  then 
averaging  over  the  fast-varying  gyrophase  a.  Here.  en  =  B/B.  .Assuming  \VKB 
ansatz  for  the  variations  perpendicular  to  B.  i.e..  <5f(x.  v.  t)  = 


6i{t, v)exp(ij"  J(k_L*dx1  —  iwt)  with  t  being  the  coordinate  along  B,  we  then  have, 
for  isotropic  plasmas  considered  here, 


ff-  (q/m)[^^--  J0(X)«5tyelLk(— )F0]  +  <5ge‘Lk  ,  (1) 

at  uj 

QF0  =  (ud/dt  +  u.)F0,  (2) 

i*F0  =■  (kX  eii/fi)’  VF0,  X==kj.p,  Lk  =  k_,  ■  e(|  X  v/f2  corresponds  to  finite 
Larmor-radius  effect  via  transformation  between  x  and  X,  and  the  field 

variables  are  56,  cdS-v/dt**  iw6An,  and  <©n.  Meanwhile.  Sg,  sa:  sfies  the 

following  linearized  gyrokinetic  equation 

where 

<5S  —  J  0(<5<j>  —  6u>)+  J0u ;d<5ii;/u/  +  Jivi<5Bn/k1c ,  (-4) 

wd=«kd-vd,  and  vd  is  the  magnetic  drift  velocity  given  by 

vd  —  en  X  [(vi/2)VlnB+  vfi/cj/fi  with  k  =■  (eii  •  V )en  being  the  curvature.  Given 
particles  are  magnetically  trapped  between  l\  and  l«,  Eq.  (3)  can  be  readily 
solved.  In  particular,  one  finds 

<5G  =  -^££g(cr)  -  — (— )QF0{cotI^  cos  1^  rb<5ScosI^t 

-  <j  m 

+  cosI^r^Ssinl^  +  -<5SsinI*e}  ,  (5) 

where  <7=sgn(v„),  l£  =  J^d^ui—  wd)/ 1  v„  | ,  rb  =  J^d£/ 1 v„  |  and 
Aa  [/^'Ad^/jv,,  |]/tv  Note  that  cot  1^  *  £  (w  —  wd  —  Ku>b )“ 1  /r0 ,  where 
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u>b » n /t^.  Equation  (5),  thus,  contains  wave- particle  interactions  via  the 
magnetic  drift  and  bounce  resonances. 


For  hydromagnetic  waves,  the  relevant  field  equations  are  the  quasi¬ 
neutrality  condition  and  the  Ampere's  law.  Specifically,  the  quasi-neutralitv 
condition  yields 


£(q2M  <  -  (“)F0<5t*/Jo  +  ( — )J 0<$G  >  ]  -  0 ,  (6) 

^  ot  jj  q 

with  j  =  species  and  <...>  =  4ir  f  (...)  Bd^de/jv,||.  The  perpendicular 
.Ampere's  law  gives 

5B„  -  <  qvaJi(— — F0J0^-  SG)  >, .  (7) 

cKj_  m  uj 

Meanwhile,  the  parallel  -Ampere’s  law  when  combined  with  the  quasi-neutrality 
condition  yields  the  following  generalized  vorticity  equation 


_1 c-  <9  kl  d 

47T  jj-  dt  B  91 


div  =* 


£q,  <  —  -MG  +  -2-[<5o 
,  jj  m 


d¥0 

de 


-  — F0(J^o  +  J‘  —  6w  +  J,  J,  T~~t5Bii )]  >, . 

'jj  jj  K  C 


(3) 


Equations  (6)  to  (8)  along  with  <5G  given  by  Eq.  (5)  form  the  complete  set  of 
equations,  which,  however,  generally  requires  numerical  solutions.  Subsidiary 
expansions  are  needed  in  order  to  make  analytical  progress  and  they  are  carried 
out  in  the  next  section. 
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EIGENMODE  EQUATIONS  AND  GENERAL  STABILITY  PROPERTIES 


Consistent  with  observations  at  the  geosynchronous  orbic,  the  ring-current 
plasma  of  electrons  and,  mainly,  protons  is  taken  to  consist  of  core  (c),  warm 
(w),  and  hot  (h)  components.  The  corresponding  temperatures  are,  formally. 
Tc  ~  0(1)  keV,  Tw  ~  0(10)  keV,  and  Th  ~  O(102)  keV.  Thus. 
Tc/Tw  -  Tw/Th  ~  10~  1  =  Sr  and  6  is  the  smallness  parameter  for  the  subsidiary 
expansions.  Meanwhile,  we  have  n~O(10/cm3,  which  takes  into  account  the 
presence  of  heavy  elements  and  B~  0(102)^.  In  terms  of  <5,  we  also  assume,  for 
recovery  phases.  d0  —  ~  <5  and  db  ~  62,  which  implies  nw/nc~62  and 

n^/h,.  —  <5°.  As  to  frequency  and  wavelength  orderings,  we  assume 
0(Va/L0),  k||L0  —  0(1),  and  m~  O(102).  We,  therefore,  have  k ±pih  ~  0(1), 
kj.Piw  ~  0(6),  kj./jjc  ~  0(6"),  and  kj.pe  being  negligible.  In  addition,  we  have 
-x’.c/w~  6o/2;  w,w/w~  6i/2;  ^dh/^~  <5~3/2;  with 

jj*  ~~  ew.  being  the  typical  diamagnetic  drift  frequency,  and  u/b  1C/ 

■^b.iw/w ~  6~ l//2,  *ub ilb/uJ—  (5~3|/2  as  well  as  -u/u; be  being  negligible.  Briefly 
stated,  the  physical  significance  of  the  above  orderings  are  (i)  while  the  core 
component  supports  the  hydromagnetic  waves,  the  warm  and  hot  component 
provide  the  instability-driving  finite-d  effects,  (ii)  while  the  core  ions  respond  to 
local  perturbations,  the  rest  of  the  components  mainly  respond  to  bounce- 
averaged  perturbations,  (iii)  for  electrons,  only  the  warm  component  experiences 
wave-particle  interactions  via  the  magnetic  drift  resonances,  (iv)  for  ions, 
however,  while  the  wave  component  experiences  both  the  magnetic  drift  and  the 
bounce  resonances,  the  hot  component  experiences  magnetic  drift-bounce 
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resonances,  and,  finally,  (v)  only  the  finite  Larmor-radius  effects  of  hot  ions 
need  to  be  kept. 


Application  of  the  above  orderings  further  simplifies  the  set  of  equations 
derived  in  the  preceding  section.  Specifically,  to  the  lowest  order  of  accuracy, 
Eq.  (6),  the  quasi-neutrality  condition,  reduces  to 

dF0  -  —  3F0  - - 

<  “a — (<5£—  <5£)>e  =  < — {ud5o—  -jd<5o)>e ,  (9) 

at  at 

where  ^j(6<p  —  <5u>)  and  is,  thus,  proportional  to  the  parallel  electric  field, 
<5E||.  Equation  (9)  indicates  that  in  realistic  magnetic  field  geometries.  <5En  could 
be  generated  via  magnetic  drifts.  Equation  (7),  the  perpendicular  Ampere’s  law. 
meanwhile,  becomes 


<©n 


4a- 


B 


wB 


2  {Qi(-lj(u;*i*.i>  '■x>*  ,te 


<9F0  v 


3F, 


+  qe  <  — - —  (ujdSa+  5£—  <5£)>e  +  q,  <  — - r-^d'5®>iw}  • 

at  2  at  2 


(10) 


where  n,u>,  t,  —  <(v}/2)  ^«F0 >, ,  njWKit|  —  <(vi/'2)wKF0>] ,  and  =* 
(ei|X/c)  ■  k/flj.  Note  that  Eq.  (10)  corresponds  to  the  magnetostatic 
approximation,  i .  e . ,  B<5Bn  +  da-dTA  =  C.  Finally,  the  generalized  vorticity 


equation  becomes 


.  r>  0,0  *  o 

3  ki  3  *  cj“  ^  Sttii,  e“ 

Ba7  T~aeH+  ~tTs° - -“..«)«» 


47re*  „  <^F0  „  f-  ^ 

- -<  “a — <*>3  W  >ic 

mjC*  c'€ 

47re2  <9F0 

+  i£V<^7*Wd^>e 

mec“ 


o  c-  4n  _  ,  q*  .  ^  a*  o  — tt  ^ 

3^>ic - -  £  (~)i  <  “oT^d^d^  >, 

c"  ]  — e,iw  m 


+  1—7-  E  (a-)j<QF0E^-3;d-Ka4)K'dJ0coSl|lJ0u;d6<5cosl|i  .  (11) 

c"  j— e.i  m  K 

Here,  we  note  that  <5Bn  does  not  appear  due  to  the  low-J  ordering,  j—  6.  and 
that  <5f  is  related  of  6<t>  via  Eq.  (9).  Equations  (9)  to  (11),  thus,  represent  a  set  of 
reduced  equations  more  susceptible  to  analytical  and  numerical  treatments. 

Some  interesting  insights  to  the  general  stability  properties  can,  however,  be 
obtained  from  the  quadratic  forms  readily  derivable  from  Eqs.  (9)  and  (11). 
Denoting  |  |  a|  |  ==  j>  ^  a  d£/B  with  ^s  and  £\j  being  the  end  points  of  the  field 

A  O  O 

line  and  5<p(£=  t$,  £y)=>  0,  we  have,  for  c*r  =2  up  +  2i  -•, ur , 


«?l|ki|6i|2/Vi| 


1  k'i  |  a<5o/a^  i  -  i  j  ~r3, 


wjWt- I*/2|  |ki|«i|2/Vi|  I; 


where 


(87rn,e2/mic2)|  |  vKX (w«  tl  -  «/.  te )  |  ho  |2  |  | 

+  (16rr2/c2)/d/id€{  £  [(q2/m)(dF0/de)rb  j ~d<5o  j 2], 


+  [(q2/m)(dF0/de)rbu4  |<5®|2l.c  +  [(qVm)(<5F0/3e)rb(  |  <5^  | 2  -  |<5f|2)|e}. 
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li-  (15) 


h**  i&Tr~Jdnde  £  El'^”TbQFo^(w—  wd— Kwb)  |wdJo60cosI* 

C"  j-e,«  K  m 

We  now  discuss  physical  meanings  of  the  various  terms.  For  typical  cases  where 
and  <9F a/dt<Q.  the  first  term  in  Rd  corresponds  to  the  bad- 
curvature  interchange/ballooning  destabilization,  while  the  second  and  third 
'erms  correspond  to  stabilization.  For  5<b  =■  constant  and  balancing  these  three 
terms,  one  can  recover  the  well-known  threshold  condition  on  the  steepness  of 
the  pressure  gradients  for  the  interchange  stability.  On  the  other  hand,  noting 
|  <5£|"—  |  <S^ | 2  >  0  from  Schwarz’  inequality,  the  last  term  in  corresponds  to 
stabilization  due  to  finite  electrostatic  coupling.  Due  to  the  present  d~0(<5)<  1 
ordering,  destabilizing  (stabilizing)  R^,  i.e.,  Rp  >Q(<0),  only  reduces  (increases) 
cjr  and  has  no  direct  bearing  on  the  stability  properties.  Extension  of  the 
present  theory  to  a?-"-0(1)  can  certainly  be  conceived  and.  in  such  cases,  the 
corresponding  Eq.  (12)  could  be  expected  to  play  the  role  of  a  generalized 
kinetic  energy  principle.  Thus,  in  the  present  work,  as  we  have  emphasized 
earlier  and  indicated  by  Eq.  (13),  stability  properties  enter  only  via  wave- 
particle  resonances  expressed  by  of  Eq.  (15).  To  go  further  with  the  stability 
analyses,  we  need  to  specify  QF0,  i.e.,  both  the  spatial  and  velocity  dependences 
of  F0.  This  is  carried  out  in  the  next  section  assuming  adiabaticaily  injected 
energetic  particles. 

STABILITY  OF  AD IAB ATIC ALL Y  INJECTED  PARTICLES 

Noting  that,  in  terms  of  the  longitudinal  invariant  J,  QF0  can  be  expressed 
as 
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(16) 


QF0  -  (w-«d)-^-| m.l  +  (kXe„/fi)  •  VL-^  |  M,  j . 

Here,  L,  the  invariant  magnetic  latitude,  is  the  equatorial  distance  of  the  field 
line  from  the  symmetry  axis.  Thus,  assuming  the  energetic  ring-current 
particles  are  adiabatically  injected  from  the  tail  such  that  <9F0/3L  |  M  j  =  0,  we 
then  have  QF0  *■  (cj— wd)<9F0/3£,  negligible  electron  contribution  Ip  e  =  0  since 
magnetic  drift  resonances  conserves  /i  and  J,  and 

=*  — ~  (“  )i/d*ide(-^)£k;r6(wr--Td-ku;b) 
c*  m  ae  ^ 

X  |cjdJ0cosI^<5|“ .  (17) 


At  the  geosynchronous  orbit,  where  proton  pressure,  typically,  decreases 
radially  outward  ( dP,/3L<0 ),  we  have  w*,wdi>0  and,  since  T„F0  —  Tdl 
dF 

<9F0,/3e,  - 2L  <  0.  Equation  (17),  thus,  indicates  that  resonant  particles 

de  . 

M.  L 

contribute  to  destabilization  (stabilization)  if  K-;r<0(>0).  For  the  K-jt  <  0 
destabilizing  effect,  the  resonance  condition,  u iT  —  Td  —  Kwb  =**  0.  becomes 


Ku»b 

wr 


1  +  u^b  |  —  |  ■ 

'Jjr 


(18) 


Equation  (18)  indicates  that  (i)  only  hot  protons  with  |Td/u;b  |  lb  —  kj^h-Si  |KJ 
contribute  to  the  destabilizing  resonances,  and  (ii)  1  >  wr/wd  lb  >  0,  i.e..  the 
excited  instabilities  have  azimuthai  phase  velocities  in  the  same  direction  but 
smaller  than  vd  lh,  which  is  westward  propagating  for  energetic  ring-current 
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protons.  To  understand  the  instability  mechanism  physically,  let  us  assume 
\u>zjZj&  |  «  1  and  the  resonance  condition  becomes  w<j  =:  |K|wb.  The  resonant 
hot  protons,  thus,  execute  radial  excursions  in  L  while  conserving  fi  and  e 
(instead  of  J).  Noting  that  dFQl/dL\  ^  £  <  0,  the  instability  is,  therefore,  fed  by 
the  hot  ion  pressure  gradient.  Noting,  furthermore.  dF0X/de  <0  and  the  j£ 
factor  in  Eq.  (17),  the  jkj^Pih  |  ^  |K|  resonance  constraint  then  implies  that  the 
K=»— 1  resonance  term  is  the  most  important  one.  For  a  dipole  B  symmetric 
with  respect  to  the  geomagnetic  equator  at  £  =»  0  (i.e.,  ty  =*  —  £?  and  £2  —  ^t)>  we 
can  easily  see  from  Eq.  (17)  that  only  odd  So  responds  to  the  K=»  —  1  resonance. 

Assuming  hydromagnetic  waves  with  odd  So  are  resonantly  destabilized  by 
the  hot  ions,  let  us  now  examine  the  stabilizing  effects,  which  require  Ku>r  >  0. 
The  resonance  condition  then  becomes  1  *  uJ<j/u;r  +  |K/wr  I  aQd  Tj/w.  >  0- 
Thus,  only  the  warm  ions  contribute  to  the  stabilization  via  the  bounce 
resonances.  Since  in  the  presence  ordering  |u//wb  lw|  ~  S1^,  the  K  =  1  term 
dominates. 

To  assess  the  relative  importance  of  the  hot-ion  destabilization  and  warm-ion 
stabilization,  we  further  assume  F0  for  both  warm  and  hot  ions  can  be 
represented  by  a  single  Maxwellian  distribution:  F01  =  n,w(2Te0)“ 
exp(— e/€0),  where  e0(L)  is  the  characteristic  energy  of  ring-current  protons. 
Assuming  kd  k0  with  X  en  in  the  westward  direction  and  keeping 

only  the  K  »»  +1  terms,  Eq.  (17)  becomes,  approximately, 


(19) 


I*  -  Ift  *  ft  l^|2(^){J3(l)(k,oPor7exp[-(k,0p0)-2] 

-  (— )7  exp[-(wr/wbo)2]}, 
wbo 

where  p0  =■  y/2to  / Qj0 ,  and  k<j0  and  B0  correspond  to  equatorial  values.  Noting 
that,  from  Eq.  (12),  u ;r  is  independent  of  k$,  maximizing  the  destabilizing  term 
in  Eq.  (19)  then  leads  to 

k*mp0  ~  0-5.  (20) 

We  remark  that,  for  ring-current  protons  at  the  geosynchronous  orbit,  we  have 
typically  e0  —  30  keV  and  Eq.  (20)  then  predicts  m  ~  0(102),  which  is 
consistent  with  satellite  observations.  With  k^  25  k^m,  Eq.  (19)  shows  that  the 
instability  threshold  condition  is  essentially  given  by 

I— I  <  1-  (21) 

Meanwhile,  from  Eq.  (13),  the  growth  rate  is  estimated  to  be 

«i  ~  /?i(V.4o/L)7wr.  (22) 

Both  Eq.  (21)  and  .Eq.  (22)  suggest  that  the  most  unstable  mode  is  an  odd-<5cz> 
eigenmode  with  the  lowest  u/r.  Now  uT  is  given  by  Eq.  (12)  along  with  Eqs.  (9) 
and  (14).  Since  |R^|  ~  0(/?)  ~  0(6)  <  1,  the  lowest-wr  odd-60  mode  tends  to 
be  the  lowest  antisymmetric  mode,  i.e.,  the  second  harmonic.  However.  Rj,  can 
make  significant  corrections  to  wr  and  it  is  worthwhile  to  evaluate  R^  for  the 
case  of  Maxwellian  adiabaticaiiy  injected  particles  and  odd-<5o  eigenmodes. 
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First,  we  have  from  Eq.  (9) 


<5£  =“  —  — <WdF0  >ec  i5o .  (23) 

Qoe 

Correspondingly,  Rg  is  then  given  by 

Rg  -  (8ff/c2)  2  [(nqVm)l  I  5<j>  | 2  |  |], 

j««ew,  iw 

+  [(47rq2/mc2€0)|  |(<i3d  -  uji)<jiF0>  |<5ci>|2  j  |]ic 
+  [(4;rq2/mc2e0)  |  |  ( <-JdoJdF0 >  —  <^dF0>  Vno)  I  <5<Z>  1 2  I  I  ]ec  •  (-4) 

Now.  for  a  dipole  magnetic  field,  one  can  readily  show  that  a>d  (£)  decreases  away 
from  the  equator  and,  hence,  <(uJd  ~  ^d)wdF0>  is  positive  (negative)  away 

A  ^ 

from  (near)  the  equator.  For  an  odd  b<t>,  \8<j> J"  is  peaked  away  from  the 
equator.  We,  therefore,  expect  |  |  <(uJd  —  wa)u?dF0>  J  | 2  |  |  >  0.  Employing, 
furthermore,  the  inequality  <u;3F0>  <F0>  >  <^dF0>2  then  demonstrates 
that  all  three  terms  in  Eq.  (24)  are  positive-definite,  that  is. 

Rg  >  0.  (25) 

Thus,  in  the  case  considered  here,  both  electrons  and  protons  will  contribute  to 
a  significant  finite-d  reduction  in  ujt,  which  is  destabilizing  both  in  terms  of 
reducing  the  bounce-resonance  stabilization  as  well  as  enhancing  the  growth 
rate. 

Finally,  let  us  consider  the  compressional  component,  <5Bn,  which  from 
Eq.  (10)  is  given  by 
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<$Bii 

B 


±L.^{[iL<LL^dF0>U  +  [7- <~-F0(ud-ujd)> 

tQ  -  — 


wB 

2 

+  qe[7~<  -^Fo^d 


lie 


Note  that  <©n  is  also  antisymmetric  as  <5©.  Both  electrons  and  core  ions, 
however,  introduce  additional  structures  along  B.  Neglecting  these 
complexities,  we  have,  approximately, 

|®ii|/|®l|  «  1.7 (27) 


where  we  have  assumed  the  ion  pressure  has  L~73  dependence  and  knL  ~  2  for 
the  second  harmonic  standing  wave.  Equation  (27)  indicates  that  the 
instabilities  tend  to  have  significant  compressions!  components. 


SUMMARY  AND  DISCUSSIONS 


In  the  preceding  sections,  we  have  presented  a  theory  on  the  excitations  of 
magnetospheric  hydromagnetic  waves  by  energetic  ring-current  particles.  Since 
the  theoretical  approach  employs  the  gyrokinetic  equations,  effects  due  to  finite 
ion  Larmor  radii,  realistic  magnetic  field  geometries  and  plasma  equilibria  are 
included.  By  assuming  a  low-d  and  isotropic  plasma,  the  present  work 
concentrates  on  excitations  via  wave-particle  resonances.  A  reduced  set  of 
eigenmode  equations  is  then  derived  by  further  assuming  the  plasma  consists  of 
core  and  energetic  components.  For  adiabatically  injected  particles,  the  present 
theory  predicts  the  following  specific  instability  features:  (1)  destabilization  via 
magnetic  drift-bounce  resonances  of  hot  ions,  T,i,  ~  O(l02)  keV;  (2)  the 
unstable  eigenmodes  are  antisymmetric  along  B  (i.e.,  odd  <5d>);  (3)  the  most 
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unstable  eigenmode  tends  to  have  k0mpo  2:  0.5  and,  correspondingly,  the  lowest 
real  frequency;  (4)  the  corresponding  azimuthal  phase  velocity  is  in  the  same 
direction  (i.e.,  westward)  but  smaller  than  the  hot-ion  magnetic  drift  velocity; 
(5)  the  instability,  in  general,  has  significant  transverse  and  compressional 
components;  and,  finally  (6)  finite-, 3  (plasma  pressure)  values  can  significantly 
reduce  the  real  frequencies.  These  features,  we  note,  are  consistent  with 
satellite  observations  mentioned  earlier. 

The  apparent  success  of  the  present  theory  has  obviously  motivated  us  to 
examine  the  more  general  case  of  a  3  —  0(1),  anisotropic  plasma,  where  one. 
perhaps,  can  hope  for  a  description  unifying  destabilizing  mechanisms  due  to 
both  velocity-space  anisotropy  (i.e.,  the  drift  mirror  mode)  and  wave-particle 
resonances.  This  attempt  for  a  unified  description  is  currently  in  progress  and 
the  results  will  be  presented  in  a  future  publication. 
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SECTION  H:  PLASMA-WALL  PHYSICS,  THEORY  AND  SIMULATION 


A.  Multi-Scale  Particle  Simulation  of  Bounded  Plasmas 

S.  E.  Parker  and  C.  K.  Birdall,  UC-Berkeley, 

A.  Friedman  and  S.  L.  Ray,  LLNL 

Our  basic  goal  is  to  study  the  kinetic  behavior  of  a  bounded  plasma,  including 
the  boundary  layer  (or  sheath)  effects.  We  want  to  understand  how  the  bulk  plasma 
behavior  with  long  space  and  time  scales  interacts  with  the  boundary  layer  where 
short  space  scales  (~  A q, ),  and  time  scales  ( ~  ~el )  are  important.  Currently, 
explicit  codes  have  been  successful  at  modeling  the  sheath  dynamics  for  system 
lengths  less  than  100  \oel-  These  codes  have  a  constraint  on  ^peAf  for  stability 
and  accuracy  which  restricts  them  to  short  systems  and  times  The  newer  Direct 
Implicit  Codes2,3  have  been  successful  at  relaxing  the  w-peA t  constraint  for  stability, 
but  are  restricted  by  accuracy  to  small  kvre^t.  where  k~l  is  the  characteristic  scale 
length  of  the  electric  field.  Typically,  fc-1  would  be  the  order  of  Ape  in  the  sheath 
region,  hence,  rewriting  kv 7\.At  =  (k\pe){^jpe^t),  we  see  that  short  time  steps  are 
necessary  to  resolve  the  sheath.  Thus,  if  one  needs  to  model  the  boundary  layer 
accurately,  the  implicit  method  has  no  advantage. 

The  constraint  on  a  individual  particle  with  velocity  v  is  bAi.  This  assures 
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proper  sampling  of  the  electric  field.  Since  the  most  severe  restriction  on  kvAt  is 
local  to  the  boundary,  we  can  take  advantage  of  this  and  move  parti.les  at  different 
time  scales  depending  on  their  location  in  phase  space  (x,v).  We  term  this  the 
Multi-Scale  Method?.  The  majority  of  the  electrons  are  in  the  bulk  plasma  region 
where  nonuniformity  in  field  and  density  are  mild  compared  to  the  sheath  region. 
All  of  the  ions  and  most  of  the  electrons  can  be  moved  much  less  frequently  than  the 
minority  of  fast  time  scale  electrons.  By  allowing  the  time  step  at  which  particles 
are  moved  to  vary,  we  can  greatly  reduce  computing  time  and  permit  modeling  of 
larger  systems. 

Our  model  is  one  dimensional,  electrostatic,  unmagnetized,  and  collisionless 
t  We  intend  to  include  Coulomb  collisions  in  the  future).  The  boundaries  are  ab¬ 
sorbing  and  there  is  a  prescribed  distributed  source.  S(x.v:t).  One  motivating  ap¬ 
plication  is  the  flow  of  a  plasma  along  magnetic  field  lines  to  a  conducting  wall. 
This  problem  is  important  for  the  understanding  of  how  plasma  flows  to  the  >li- 
vertor  plate  or  limiter  in  a  Tokamak  fusion  device0.  One  is  interested  in  how  the 
boundary  layer  effects  the  bulk  plasma. 

Current  Progress 

We  have  developed  a  bounded  one  dimensional  electrostatic  particle  simulation 
that  uses  the  Direct  Implicit  Method6.  A  consitant  way  of  prescribing  the  boundary 
conditions  when  using  the  Direct  Implicit  Method  has  been  developed.  Test  runs 
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show  that  the  sheath  is  not  properly  resolved  which  motivates  implementation  of 
the  Multi-Scale  Method. 

A  periodic  code  called  MIST  (Multi-Scale  Implicit  Testbed)  is  being  developed 
to  test  the  Multi-Scale  Technique’.  The  algorithm  has  been  worked  out  and  is 
explained  in  Reference  4.  One  of  the  current  research  issues  is  the  criterion  for 
choosing  a  particles  At,  see  Reference  7. 
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B.  Numerical  Error  in  Electron  Orbits  with  Large  ^ceA t 

S.  E.  Parker  and  C.  K.  Birdsall 
Electronics  Research  Laboratoy,  U.C.  Berkeley 
Berkeley,  California  94720 

Inw  duction 

We  have  been  using  a  bounded  2d  electrostatic  particle  simulation1  to  study 
electrostatic  effects  induced  by  nonuniform  magnetic  fields*.  It  was  found  that 
running  with  large  values  of  -jceA/  did  cot  se»m  op  e  any  significant  problems 
in  our  model.  The  wot1.:  presented  here  is  to  acldi.  ss  why  p'.asma  simulation  wuth 
large  ~jceAt  under  controlled  circumstances  can  give  good  resrdts. 

When  simulating  a  magnetized  plasma  with  electron  inertial  effects  included, 
the  limiting  time  scale  in  many  situations  is  the  electron  cyclotron  period,  with  a 
constraint  on  ,jceAf.  If  one  is  not  interested  in  short  time  scales,  this  restraint  can 
make  simulation  of  other  effects  at  longer  time  scales,  e.g.  ion  transit  times,  very- 
expensive.  Particle  simulation  has  been  done  using  the  guiding  center  equations  for 
electrons3  to  relax  the  *oce  A<  constraint.  Alternately,  we  are  exploring  the  feasibility 
of  running  with  large  ^'C(.Af  using  the  standard  integration  scheme  of  Boris4.  Wp 
begin  by  comparing  single  particle  orbits  with  large  u>ceA t  to  those  with  small 
nceA t.  Then,  the  error  introduced  when  nC(.Af  is  large  is  discussed. 


Results  from  a  single  particle  mover  for  uiceAt  1 


To  study  the  errors  introduced  with  large  -nce  At,  we  used  a  single  particle  mover 
with  two  spatial  coordinates,  (x,y)  and,  three  velocity  coordinates,  (u£ ,  vv ,  vz ).  The 
following  parameters  are  the  same  in  all  the  runs:  v(t  =  0)  =  (0.1,  0.0, 0.4)  and. 
■i.  =  1.  The  magnetic  field,  electric  field,  and  time  step  is  varied.  Figure  1  has 
two  runs  with  ajce  A t  =  0.5  and  50;  B  =  (250,  0,  0),  and  E  =  (0, 0, 1)  for  both  runs. 
V£)£  =  0.004  y  as  predicted.  For  the  run  with  ,nceAt  =  50;  Tg,eff  ~  t'j_A t.  and 
wc  eff  ~  0.04  as  expected  from  the  discussion  below.  Figure  2  has  twro  runs  with 
a,-ceA t  ^  0.5  and  50.  In  these  runs  the  magnetic  field  is  nonumform  and  there  is  no 
electric  field.  Note  that  the  averaged  orbits  are  approximately  equal. 

Single  particle  orbits  with  ^ceAt>  1 


In  this  section  we  characterize  the  electron  orbit  when  -*ce  At  is  large.  The 
v  x  B  integration  scheme  gives  the  angle  of  rotation  for  one  time  step  as4: 


9  =  2  arc  tan 


(1) 


As  ,nceAt  -*  x,  9  — *  it.  The  orbit  for  ,*;ceAt  >■  1  can  be  described  as  rapid  bounce 
back  and  forth  motion  with  the  perpendicular  velocity,  and  a  slow  precession  about 
the  guiding  center  axis.  For  large  ,nceA t,  the  angle  of  rotation  is  approximately: 


9 


7T  — 


u^re  At 


(2) 
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The  precession  is  the  result  of  the  rotation  angle  being  a  small  amount  8  less  than 
7T,  6  «  n  —  9.  The  frequency  of  this  precession  is: 


At  wceAf 


The  “effective”  cyclotron  frequency  is  just  or: 


1  (  4 

-  7T - 

At  \  uJc,At 


(3) 


(4) 


for  large  ~uceAt.  The  "effective”  gyroradius  for  all  values  of  -*ceA t  is  given  by: 

(51 

Hence  ~  O(At),  and  can  be  made  small  if  it  is  a  major  cause  of  noise  in  the 

simulation  model.  Note  that  u jc  tff  goes  to  ~  for  large  -;ceAt. 

Guiding  center  motion  with  o,ceAt>  1 

Let  us  analyze  the  guiding  center  motion  to  see  what  error  is  introduced  when 
^ceA t  is  made  large.  We  start  with  the  integration  scheme  eliminating  v  using 
v"-2  =  (xn  -  xn-1)/Af  and,  vn+i  =  (x"+1  -  xn)/A t: 

-^■(xr,+1  -2xn+xr-1)  =  ^-{E(x")  +  ^(x"+1  -  xn_1)  x  B(xn)|  (6) 


Vj_At 
!  sin  ( !■£) 


1  -U 


(^)! 


Let  the  magnetic  field  B  be  in  the  x  -  y  plane  and  define  the  following  coordinates: 
ei  =  §.  e2  =  z.  and  e3  =  ej  x  e2.  Let  x"  =  x£  -f  x".  where  x"  =  rg -e/ f{sin®e2  + 
cosoes),  and  <z>  =  <z>o  +  At  u>c  tff(tn).  xq  is  the  slowly  varying  guiding  center 
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motion,  and  Xi  is  the  fast  gyration.  Assume  that  ^c.e/ /  does  not  van-  much  between 
time  steps  so  we  can  use  the  approximation  %  ^c,«//(fn+1  )•  Taking  a 

Taylor  series  expansion  of  equation  (5)  ahout  Xo,  and  then  equating  the  slowly 
vaxving  terms  gives  the  guiding  center  equation: 


--L(xr-2x;  +  x;-,)  =  X{EW)  +  -l.( 


,n-rl 


m 


2A  t 


—  x 


n  —  1 


)  x  B(x” )  j  —  - 


Ueff 


m 


vb(X; 


where  Me//  =  -jfn(~‘A,e//-kO-  Using  equations  (3)  and  (5).  we  get  m,//  = 

=  M-  The  guiding  center  equation  of  motion  given  by  Northrup^  is: 


d-x  o 
dt2 


x  B(x0)j  -  -VB(X o) 
at  1  m 


(S) 


Equation  (7)  is  just  a  central  difference  approximation  to  equation  (S).  Hence,  we 
expect  the  same  first  order  drift  motion. 

Conclusions  and  future  work 


We  conclude  that  even  with  1  the  Boris  algorithm  still  gives  the 

correct  guiding  center  motion  for  a  single  particle.  This  paper  did  not  take  into 
account  the  collective  error  introduced  by  a  large  number  of  particles.  We  need  to 
examine  collective  numerical  instability. 

A  possible  improvement  that  would  eliminate  the  and  r}eff  effects  on 

the  fields  would  be  to  subtract  off  the  rapid  bounce  motion  before  the  wieghting  of 
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particles  to  calculate  the  charge  density. 


(9) 


Where  the  /  subscript  represents  the  position  for  particle  wieghting  for  the  charge 
density.  It  would  only  apply  when  cjceAi  1.  This  is  just  a  suggestion  and  needs 
to  be  tested  as  part  of  our  future  work. 
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C.  Electron-Neutral  Particle  Collisions  in  ESI 


Julia  Little 


1.  Introduction 

We  axe  adapting  ESI  to  incorporate  electron-neutral  particle  collisions.  In  the  for¬ 
mulation  of  the  problem,  we  assert  that  the  collisions  are  elastic  since  the  mass  of  the 
electron  is  negligible  when  compared  to  the  mass  of  the  neutral  particle.  We  note  that 
the  radiation  process,  considered  to  be  negligible,  is  omitted1.  We  also  note  that  electron- 
molecule  collisions  are  more  complicated  than  electron-atom  collisions1  ,  but  we  treat  only 
the  electron-atom  collisions. 

The  electron-neutral  collisions  are  modelled  with  electrons  scattering  off  a  fixed  back¬ 
ground  of  neutrals  as  if  the  neutrals  are  hard  and  massive  spheres.  Scattering  angles  are 
chosen  uniformly  randomly  in  the  interval  0  to  t  and  phase  angles  are  chosen  uniformly 
randomly  in  the  interval  0  to  2x.  The  details  of  the  simulation  of  hard  sphere  scattering 
and  of  the  method  of  deciding  when  an  electron  will  scatter  are  described  in  the  next  two 
sections.  In  the  fourth  section,  we  discuss  how  the  computer  routine  modelling  electron- 
neutral  collisions  in  the  above  manner  was  tested  and  showed  some  expected  statistics  for 
the  angular  spread  in  time  of  a  beam  of  virtually  ballistic  particles.  Finally,  we  will  suggest 
further  improvements  and  uses  of  ESI  with  the  elastic  collisions. 

The  model  above  does  not  handle  charge  exchange  between  species  or  ionization  and 
excitation  of  neutrals.  However,  we  believe  it  is  a  good  first  step  for  the  effect  of  elastic 
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collisions  on  tlie  plasma  behavior  in  the  regime  where  electron- neutral  collisions  exceed 
Coulomb  collisions,  as  is  true  for  weakly  ionized  plasmas. 

2.  The  Model  and  Simulation  of  Hard  Sphere  Scattering 

We  essentially  use  the  scattering  routine  written  in  1985  by  Perry  Gray  for  use  in  PDWl 
but  with  minor  changes.  Since  ESI  is  a  one  dimensional  program,  only  the  x  component  of 
velocity,  vx  is  used  in  finding  the  new  position  for  a  particle.  (If  the  plasma  is  magnitized. 
then  the  y  component  of  the  velocity,  vy  is  needed  for  the  v  x  B  force,  but  again  only  vx  is 
used  to  advance  the  particle’s  position.) 

In  our  program,  however,  when  an  electron  is  to  collide  with  a  neutral,  we  include  vv 
and  introduce  a  z  component  of  velocity.  vz.  (See  Fig.l)  (Initially,  we  set  cy  to  0  for  all 
particles  if  the  plasma  is  unmagnitized.  Otherwise,  the  determination  of  the  initial  values 
of  vy  is  left  unchanged.  We  give  vz  a  Maxwellian  distribution  if  the  plasma  is  warm  and 
otherwise  set  vz  to  0.)  Prior  to  making  a  collision.  ESI  assumes  the  velocity  vector,  v  has 
components  on  a  sphere  of  radius  |u|,  where  |vj  is  the  magnitude  of  v: 

vx  =  |u|cosa, 

vy  =  |v|sinacosd,  and 

vz  =  |u|sinasind. 

Then  the  scattering  angle,  7  and  the  phase  angle,  6  axe  both  chosen  uniformly  randomly 
in  the  interval  0  to  jt  and  0  to  27T,  respectively.  Hence  after  the  collision  the  new  values  for 
the  velocity  components  become 

=  Mc°s7, 
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vycoU  =  |y|sin7cos^,  and 


Vt'oit  =  Msimrsintf. 

Note  that  the  magnitude  of  the  velocity  of  a  particle  remains  constant  for  all  time  steps. 
Hence,  the  effect  of  an  electron-neutral  collision  is  the  rotation  of  v  through  a  uniformly 
random  angle  about  a  sphere  of  radius  |w|.  Only  the  direction  of  v  changes  and  not  its 
magnitude. 

3.  How  The  Decision  Of  Which  Electron  Is  To  Collide  Is  Made 

Given  a  physical  system  of  electrons  and  neutrals  including  electron-neutral  collisions, 
an  upper  estimate  of  the  collision  frequency,  v  must  be  made  in  advance  in  order  to  appro¬ 
priately  determine  the  time  step.  At  to  be  used  in  the  simulation.  For  example,  one  may 
want  approximately  all  particles  to  collide  in  100  time  steps.  Given  v  by  the  physics  of  the 
system,  then  At  is  determined  from  the  equation  i/A t  =  (Actually,  one  approximates 
the  relative  velocity  of  the  electron  with  the  neutral,  ure(  and  the  mean  free  path  for  the 
svstem.Am/„  and  solves  At  =  rka  for  At,  since  v  =  ,<'Mi  .  Consequently,  a  convenient 

Am/p  100  *m/ p 

input  parameter  for  ESI  with  elastic  collisions  is  the  inverse  of  the  mean  free  path.  - l--  . ) 

Having  chosen  At,  we  are  now  able  to  assign  a  probability  of  collision  to  an  electron  in 
such  a  way  that  the  probability  of  the  electron  colliding  is  a  function  of  its  velocity.  If  the 
magnitude  of  its  velocity  vector,  ju|  is  large,  the  probability  of  collision  will  be  high  and  if 
the  magnitude  is  Small,  the  probability  of  collision  will  be  low.  That  is.  for  each  particle 
and  at  each  time  step  we  choose  a  uniformly  random  number,  r  and  if  r  <  -r-1 — |u|Af.  the 

Am  f  p 
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particie  will  collide  according  to  the  scattering  procedure  described  in  the  previous  section. 
Otherwise,  the  particle  will  not  collide  and  no  adjustments  will  be  made  to  vx ,  vy,  and  vz. 
4.  Testing  The  Collision  Model  On  a  Beam  Of  Electrons 


To  test  ESI  with  electron-neutral  colisions,  we  made  4  runs,  each  of  which  started  out 
with  a  mono-energetic  beam  of  electrons  (at  time  zero  vx  =  1,  vy  =  0  and  oz  =  0).  In  all 
runs  the  plasma  frequency,  -jp  was  set  to  10 -12  so  that  the  measure  of  any  fields  would  be 
negligible.  The  effect  was  to  make  the  electrons  virtually  chargeless  and  hence  to  make  the 
electron-neutral  collisions  ballistic.  Also,  for  all  runs  the  length  of  the  system.  I  was  1.  The 
four  runs  differed  by  the  following  parameters: 


We  feel  confident  that  the  scattering  routine  works  effectively  in  ESI.  The  routine 
remains  to  be  tested  on  a  plasma  with  fields  and  also  for  large  uAt  or  where  l  <<  \m;p. 


Also,  the  procedure  for  deciding  whether  or  not  an  electron  with  a  given  velocity  collides 
eventually  can  be  made  vectorizable. 
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X 


Figure  1.  Velocity  angles  for  electron-neutral  elastic  scattering.  The  initial  velocity 
is  v  and  the  frnal  velocity  'is  v:ol! ,  same  magnitude  by  different  angle. 
The  scattering  is  "larg?  angle,"  not  "cumulative  small  angle",  as  with 
Coulomb  collisions. 
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d.  Particle-In- Cell  Simulations  of  Radiofrequency 
Heating  (ECRH)  of  a  Simple  Mirror  Plasma 


R.  J.  Procassini,  J.  C.  Cummings 
and 

B.  I.  Cohen  ( Lawrence  Livermore  National  Laboratory) 

Introduction 

Particle  confinement  in  mirror  plasmas  is  based  upon  the  conservation  of  kinetic 
energy  and  the  adiabaticity  of  the  magnetic  moment,  such  that  the  particles  bounce 
back  and  forth  in  the  magnetic  potential  well.  Only  a  certain  fraction  of  the  particles 
axe  trapped  or  confined  by  this  potential  well.  The  trapped  region  in  Uj.  vs-  y'|  space 
is  defined  by  the  mirror  magnetic  field  ratio  and  the  plasma  ambipolar  electrostatic 
potential.  The  main  loss  mechanism  for  particles  in  both  simple  mirror  plasmas  and 
in  the  end  plugs  of  tandem  mirror  plasmas  is  collisional  detrapping  of  the  particles 
from  this  trapped  region  into  the  loss  cone.  Once  in  the  loss  cone,  the  particles  axe 
free  to  transit  out  the  ends  of  the  device,  unhindered  by  the  magnetic  field. 

Ions  in  the  center  cell  region  of  a  tandem  mirror  device  axe  electrostatically 
confined  by  large  positive  ambipolar  potentials  which  arise  in  the  end  plugs  of  the 
device.  In  the  simplest  model,  this  potential  increases  linearly  with  the  electron 
temperature  in  the  end  plugs. 


-38- 


It  is  possible  to  improve  particle  confinement  in  mirror  devices  through  the  in¬ 
jection  of  radiofrequency  waves  at  the  fundamental  or  some  harmonic  of  the  electron 


cyclotron  freqeuency.  This  process  is  termed  Electron  Cyclotron  Resonance  Heating 
(ECRH).  The  increase  in  the  confinement  of  center  cell  ions  in  a  tandem  mirror  is 
due  to  the  heating  of  the  plug  electrons,  which  results  in  an  increase  in  the  confin¬ 
ing  electrostatic  potential  in  the  plugs.  Of  secondary  interest  is  the  enhancement 
of  electron  confinement  in  a  simple  mirror  or  tandem  mirror  plug,  which  is  due  to 
the  preferential  increase  of  the  perpendicular  velocity  relative  to  the  parallel  veloc¬ 
ity.  This  leads  to  hot  electrons  which  axe  placed  farther  into  the  trapped  region  of 
velocity  space,  thereby  requiring  a  longer  time  to  collisionally  migrate  to  the  loss 
boundary  in  velocity  space.  This  assumes  that  the  hot  electrons  remain  microstable. 


Methodology 

The  TESS  (Tandem  Experiment  Simulation  Studies)  code  is  being  developed  to 
investigate  thermal-barrier  formation  and  particle  confinement  in  a  tandem  mir¬ 
ror  plasma.  TESS  is  a  1D-3V  (one  spatial  dimension,  three  velocity  component) 
guiding-center,  relativistic,  electrostatic  particle  code,  which  includes  a  Monte  Carlo 
binary  particle  collision  package,  a  neutral-beam  heating  package  and  a  radiofre¬ 
quency  (RF)  heating  package,  which  is  the  focus  of  this  paper. 

The  RF-heating  package  is  based  upon  the  quasilinear  model  of  Berstein  and 


Baxter1,  which  treats  the  heating  of  mirror-trapped  particles.  The  model  is  valid 
for  mildly  relativistic  electrons  in  the  presence  of  a  cyclotron  resonant  applied  elec¬ 
tromagnetic  field.  The  model  assumes  that  the  cyclotron  frequency  is  much  higher 
than  the  bounce  frequency  in  the  magnetic  potential  well,  and  then  averages  over 
the  particle  trajectories  as  they  bounce  in  the  well.  One  of  the  main  results  of 
Reference  1  is  the  derivation  of  a  quasilinear  diffusion  equation  which  describes  the 
slow  evolution  of  the  heated  electron  distribution  function  fa  in  e  (total  energy)  .  n 
(magnetic  moment)  space. 


The  RF-heating  model  was  incorporated  into  the  particle  code  in  the  following 
way.  The  region  over  which  the  RF  power  is  being  injected  is  defined  by  a  center 
location  ZORf.  and  a  width  AZ/yr.  The  intensity  of  the  injected  RF  power  varies 
within  the  region  as 


Erf(z)  = 


ERF(ZORf.) 


■  [l  -I-  cos  ( 


2n(z  -  Za  ) 


AZrf 


^)] 


(1) 


where  Erf{Zorf  )  is  the  maximum  amplitude  of  the  RF  wave  within  the  region.  Par¬ 
ticles  which  lie  within  this  region  are  tested  to  determine  if  they  meet  the  resonance 
condition  R  =  0,  where 


tfl,  fc||P|| 

It  =  URF - 

7  7m, 


(-) 


with  uirjt  the  frequency  of  the  injected  wave,  l  the  cyclotron  harmonic  number.  f2,  = 


q3B/m,c  the  cyclotron  frequency  for  species  s,  7  =  1/  \J  l  —  v2/c2  the  relativistic 
factor,  fc||  the  wave  vector  parallel  to  the  magnetic  field,  and  pq  =  7m,ujj  the  particle 
momentum  parallel  to  the  magnetic  field.  Finally,  m ,  and  q ,  axe  the  species  rest 


mass  and  charge  respectively,  and  c  is  the  speed  of  light.  If  the  value  of  R  changes 
sign  in  one  timestep,  the  particle  has  gone  through  resonance,  and  is  therefore  given 
a  ‘‘kick”  in  energy. 


The  change  in  energy  that  the  resonant  particle  receives  is  given  by 


As  =  (^i)  At  +  ( 6DctAt)l/2RN  (3) 

where  e  =  7 m,c2  4-  qs$  is  the  total  particle  energy  an  for  electrostatic  potential 
$.  D(C  is  the  energy  diffusion  cofficient  from  the  quasilinear  diffusion  equation.  A t 
is  the  timestep  and  J?,v  is  a  uniformly  distributed  random  number  in  the  range 
[—1.1].  The  energy  diffusion  coefficient  is  a  specified  function  of  the  particle  rest 
mass,  charge,  parallel  velocity,  magnetic  moment,  gyroradius,  cyclotron  frequency 
and  the  parallel  and  perpendicular  wave  vectors1.  The  derivative  of  the  energy 
diffusion  coefficient  with  respect  to  total  energy  may  also  be  written  as 


dD"  dDu  dD" 
de  ~  de  +  “  dp 

where  a  =  (1  —  k\\v\\lunF)l B  a-nd  the  magnetic  moment  n  =  ~f2m,v]_/2B. 


(4) 


The  first  term  on  the  right  side  of  (3)  represents  the  change  in  total  energy 
due  to  dynamic  friction  or  drag,  while  the  second  term  represents  the  change  due 
to  diffusion  in  e,  fx  space.  These  are  the  same  two  processes  that  are  described  by 
the  Fokker-Planck  collision  operator.  Note  that  the  energy  increment  Ae  may  be 
either  positive  or  negative,  such  that  the  particle  gains  energy  from  the  wave,  or 
vice  versa,  depending  upon  the  direction  of  particle  gyration. 


Once  the  total  energy  increment  is  known,  we  calculate  the  change  in  the  particle 


parallel  momentum  and  magnetic  moment  by 


Ap|]  =  7  m,c2 


Vll 

uJRF  P\\C2 


(5) 


**  =  i(l  -  =  Icle.  (6) 

a  ijjRp  B 

This  implementation  of  the  model  was  used  earlier  by  Rognlien  in  the  Monte  Carlo 
RF-heated  mirror  code  MCPAT2,3. 


Simulation  Results 

In  addition  to  the  data  required  to  run  the  basic  electrostatic  particle  code,  TESS 
also  requires  the  following  data  for  the  RF  heating  model:  the  center  and  width  of 
the  RF  heated  region  ( ZORF  and  AZ/yr),  the  wave  frequency  (u >rf),  the  cyclotron 
harmonic  number  of  the  wave  (i),  the  peak  RF  wave  amplitude  ( Erf{Z0rf ))  and 
the  parallel  and  perpendicular  wave  vectors  (fc||  and  k±). 

Results  of  the  TESS  code  are  compared  to  those  from  MCPAT  below.  For  these 
runs,  the  magnetic  field  profiles  were  specified  in  both  codes,  and  the  self-consistent 
calculation  of  the  ambipolar  electrostatic  potential  was  disabled  in  TESS.  in  an  effort 
to  more  closely  compare  with  MCPAT.  which  does  not  have  this  capability.  The  runs 
are  therefore  simulations  of  a  simple  mirror  geometry,  with  ECRH  applied  at  the 
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bottom  of  the  magnetic  well.  The  RF  waves  are  injected  essentially  perpendicular  to 
the  magnetic  field  at  the  bottom  of  the  magnteic  well,  with  a  small  component  along 
the  magnetic  field.  The  input  data  for  the  runs  were  ZORF  =  250cm,  A Zrf  =  50cm. 
Ufljr  =  8.8  x  10lorad/sec,  Ehf{Zorf)  =  3V/cm,  t  =  1  or  fundamental  heating,  and 
kx  =  1.0cm-1  and  fcjj  =  0.1cm-1. 

The  electron  velocity-space  plots  for  the  RF  heated  region  from  TESS  and  MC¬ 
PAT  are  shown  in  Figure  1  below.  Note  the  increase  in  particle  density  in  the 
regions  of  high  fj.  and  low  V||,  which  is  the  indication  of  preferential  heating  in  the 
A-direction.  relative  to  the  ||-direction.  The  initial  thermal  velocity  in  each  direc¬ 
tion  is  vth  =  4.19  x  108  cm/sec,  such  that  the  most  energetic  particles  have  been 
heated  to  a  maximum  velocity  vXrna,  —  6vth- 

The  kinetic  energy’  profiles  along  the  axis  of  the  simple  mirror  are  shown  in  Figure 
2.  Note  that  the  TESS  result  is  an  instantaneous  snapshot,  while  the  MCPAT  result 
is  time-averaged.  The  MCPAT  profile  from  —150cm  to  150cm  corresponds  to  the 
region  from  100  to  400cm  (the  peaks  of  the  magnetic  field)  on  the  TESS  profile. 
Note  that  in  each  case,  the  parallel  kinetic  energy  is,  on  average,  half  that  of  the 
perpendic’ilar  kinetic  energy.  This  is  because  the  perpendicular  direction  comprises 
two  degrees-of-freedom,  as  compared  with  one  for  the  parallel  direction.  The  heating 
is  manifest  almost  entirely  in  the  increase  of  the  perpendicular  kinetic  energy’  in 
the  25cm  to  either  side  of  the  injection  center  point.  TESS  gives  a  maximum 
perpendicular  kinetic  energy  of  215eV,  while  MCPAT  yields  a  200eV  maximum 


-43- 


-44- 


a)  TESS  velocity  space  density  plot. 


Vx  ( cm/sc  c) 

b)  MCPAT  velocity  space  density  plot. 


Figure  1.  Electron  V||  velocity  space  density  plots  for  the  region  near  the  bottom 
of  the  magnetic  well,  showing  the  preferential  increase  in  perpendicular  velocity  due 
to  heating  by  applied  electromagnetic  waves. 


a)  Instantaneous  kinetic  energy  profile  from  TESS. 


T  versus  l 


Ax/a/  Position 


b)  Time-averaged  kinetic  energy  profile  from  MCPAT. 


Figure  2.  Electron  kinetic  energy  profiles  showing  the  increase  in  perpendicular 
kinetic  energy  in  the  RF-heated  region.  The  RF-heated  region  is  25cm  on  either 
side  of  the  magnetic  well  minimum. 
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value.  The  large  values  of  the  perpendicular  kinetic  energy  near  the  peaks  of  the 
magnetic  field  in  the  TESS  results  are  due  to  particles  which  are  near  their  turning 
points,  with  large  values  of  v±  and  nearly  zero  values  of  i>||.  This  effect  is  not  seen 
in  the  MCPAT  results  due  to  the  time- averaging  of  the  profile. 

The  results  from  TESS  compare  favorably  to  those  of  MCPAT.  The  next  step  in 
this  research  with  TESS  is  to  calculate  the  electrostatic  potential  self-consistentlv. 
allow  the  particles  to  collide  in  addition  to  their  diffusion  due  to  the  ECRH.  and 
thereby  determine  the  effect  that  RF-heating  has  on  the  overall  particle  confinement 
time.  An  additional  area  of  research  will  entail  a  comparison  of  the  results  from 
TESS  and  a  bounce-averaged  Fokker-Planck  code  for  the  simulation  of  particle 
confinement  in  an  RF-heated  simple  mirror.  This  comparison  will  be  based  upon 
velocity  distribution  function  resolution,  computational  cost  effectivness  and  ease 
of  operation. 
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-48- 


Vortex  Formation  and  Transport  to  the  Wall  in  a  Crossed-Field  Sheath*. 

K.  THEILHABER,  ERL,  University  of  California,  Berkeley,  CA  94720. 


Panicle  simulations  of  a  transversely  magnetized  sheath  have  been  conducted  with  the  aim  of  modeling 
plasma  behavior  in  the  vicinity  of  the  limiters  and  walls  of  a  fusion  device.  In  the  model  the  magnetic 
field  is  parallel  to  the  wall.  The  two-dimensional,  electrostatic,  bounded  particle  simulation  code  ES2 
has  been  used  as  a  tool  for  the  investigation  of  these  edge  effects.  The  simulations  show  that  the 
bounded  plasma  is  selfconsistently  subject  to  a  Kelvin-Helmholtz  instability,  which  is  driven  by  the 
non-uniform  electric  field  created  by  the  preferential  charging  of  the  boundries  by  ions.  This  instability 
is  seen  to  evolve  into  many  vortices  which  coalesce  into  large  vortices  (ey/T ,=1)  which  drift  parallel  to 
the  walls;  these  vortices  convect  particles  to  the  walls  at  an  anomalous  rate  much  greater  that  that  in¬ 
duced  by  collisional  diffusion.  Volume  ionization  of  neutrals  has  been  modeled  by  electron-ion  pair 
creation;  this  results  in  a  steady-state,  in  which  the  linear  edge  instability,  the  nonlinear  fluid  dynamics 
of  the  vortices,  and  the  nonlinear  dynamics  of  the  particles  scattered  to  the  walls  by  the  vortices  all  bal¬ 
ance  each  other.  This  steady-state  but  non-equilibrium  configuration  can  be  likened  to  Rayleigh-Benard 
convection.  Results  will  be  presented  for  the  scaling  of  transport  as  a  function  of  the  edge  parameters, 
and  for  the  characteristics  of  the  spectrum  of  the  ExB  turbulence.  A  computer-generated  movie  wil  be 
shown,  illustrating  vortex  formation  and  coalescence  in  the  sheath,  and  the  resulting  depletion  of  pani¬ 
cles. 


•Work  done  in  collaboration  with  C.K.  Birdsall.  Research  supported  by  U.S.  Department  of  Energy 
Contract  No.  DE-FG03-86ER53220  and  by  Office  of  Naval  Research  Contract  No.  N00014-85-K-0809. 


Ion  Acceleration  in  a  Plasma  Source  Sheath.  C.K. 
BrRDSALLTTT.5.  THEILIIABER,  S 7' KUHN,'  ERL,  Univ  of  Calif  . , 
Berkeley  CA  94720  -  A  planar  plasma  source  emitting  warm 
electrons  and  cool  ions  may  create  a  monotonically  de¬ 
creasing  potential  which  repels  some  of  the  electrons  and 
accelerates  all  of  the  ions,  in  a  model  which  is  the  dual 
of  the  Bohm  collector  sheath.  For  cold  ions  and  Boltz¬ 
mann  electrons,  0"=O(n.=n  )  gives  a  maximum  ion  flux  for 
a  source  sheath  potential  crop  of  e0/kT  =-0.5,  hence  an 
ion  velocity  of  sougd  speed.  Next,  invoking  a  real 
source  field,  (E(0)  >0) ,  leads  to  a  somewhat  smaller  max¬ 
imum  ion  flux  at  e0/kT  =-1.256,  with  ion  speed  of  1 .585v 
Last,  adding  an  initial  ion  velocity,  we  find  analytic  S 
expressions  for  ion  flux  and  plasma  potential  (measured 
at  0"=O)  which  agree  well  with  the  numerical  results  of 
Kuhn  (1)  who  used  the  proper  cut-off  Gaussians  for  both 
electron  and  ion  velocity  distributions.  Thus,  in  a  col¬ 
lisionless  plasma,  the  Bohm  "pre-collector-sheath"  can 
be  replaced  by  a  planar  source  field  (or  by  a  distributed 
source)  eliminating  the  need  for  postulating  a  pre-sheath 
acceleration  by  a  small  force  acting  over  a  long  distance 

(1)S.  Kuhn,  Plasma  Physics,  VoL23,p881, 1981;see  Fig. 2a. 
Work  supported  in  part  by  DOE  Contract  DE-FG03-86ER53220 
and  ONR  Contract  N00014-85-K-0809 . 


Prospects  for  Multi-Scale  Particle-in-Cell  Simulation 
of  Plasmas.*  A.  FRIEDMAN  and  S.L.  RAY,  LLNL,  and  S.  PARKER 
and  C.K.  BERDSALL,  U.C.  Berkeley  —  A  long-standing  goal  in  plasma 
simulation  has  been  a  method  which  could  treat  both  detailed  kinetic 
physics  and  smooth  large-scale  physics  in  an  efficient  and  natural  way. 
We  describe  a  proposed  technique1  which  would  be  suitable  for  strongly 
inhomogeneous  problems  involving  a  wide  range  of  space  and  time  scales, 
and  the  beginnings  of  our  investigation  into  its  feasibility.  For  such 
problems  to  be  tractable,  the  particles  in  each  part  of  phase  space  must 
be  advanced  on  their  own  natural  timescales,  which  are  set  by  accuracy 
criteria  ( e.g limits  on  wUapAt  or  kv&t).  In  the  past,  this  has  been 
deemed  impractical  due  to  the  requirement  of  a  self-consistent  field.  We 
propose  to  overcome  this  difficulty  by  combining  implicit  PIC  methods 
with  a  nonuniform  mesh  and  suitable  interpolations  in  time.  Since  most 
particles  are  not  processed  on  any  given  step,  on  suitable  problems  we  can 
hope  for  orders-of-magnitude  speed  improvements  over  explicit  codes. 

‘Work  performed  for  USDoE  under  contract  W-7405- ENG-48. 

1  A.  Friedman,  2nd  US- Japan  Workshop  on  Advanced  Plasma  Modeling, 
Nagoya,  March  1987. 
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Investigations  of  the  Pierce  Diode  Strane^  '  ttractor.” 
Wm.S.  Lawson,  E.R.L.,  Univ.  of  Calif.,  Berkeley  —  T1  o.nge  attrac¬ 
tor  discovered  by  Godfrey1  in  the  Pierce  Diode  system  near  :  =  2.85x  is 
studied.  While  the  results  differ  slightly  from  those  of  Godfrey  (his  time 
step  was  not  small  enough  for  accurate  convergence),  the  existence  of  the 
attractor  is  verified,  and  some  other  properties  are  noted.  The  changes 
in  the  attractor  upon  the  introduction  of  a  capacitor  into  the  external 
circuit  are  also  investigated.  The  results  show  that  even  a  high  capaci¬ 
tance  (1000  times  the  vacuum  capacitance  of  the  diode  itself)  results  in 
significant  changes  in  the  morphology  of  the  attractor  as  a  function  of 
a.  Further  reduction  in  the  external  capacitance  shifts  the  position  of 
the  attractor  slightly,  and  the  amplitude  of  the  electric  field  excursions 
doubles.  By  the  time  the  external  capacitance  reaches  five  times  the  in¬ 
ternal  capacitance,  the  strange  attractor  has  vanished. 

*  Work  was  supported  by  DOE  contract  DE-FG08-86ER53220. 

1  B.B.  Godfrey,  Phys.  Fluids  30(1987),  1553 


Direct  Implicit  Particle  Simulation  of  a  Bounded 
Plasma  System.*  S.E.  PARKER  and  C.K.  BIRDSALL,  U.C.  Berkeley 
and  A.  FRIEDMAN  and  S.L.  RAY,  LLNL  —  We  are  developing  a  di¬ 
rect  implicit  electrostatic  particle  code  which  models  a  one  dimensional 
bounded  system.  The  system  has  injection  and  absorption  of  plasma  at 
the  boundaries,  with  the  left  plane  grounded  and  the  right  plane  floating 
(open  circuit).  One  objective  is  to  use  large  mi/me  and  long  time  steps 
m  order  to  obtain  low  frequency  ion  behavior  in  detail,  including  that 
at  the  source  and  collector  sheaths.  While  the  implicit  method  resolves 
long  wavelengths  well,  there  are  questions  about  resolving  the  shorter 
wavelength  behavior  in  the  sheath  regions.  Hence  comparison  with  re¬ 
sults  from  an  explicit  one  dimensional  bounded  code1  will  be  presented. 
A  long  term  goal  will  be  to  use  this  code  as  a  test  bed  for  particle  simu¬ 
lation  with  variable  space  and  time  scales.  A  scheme  has  been  proposed 
where  particles  will  be  moved  at  their  respective  time  scales  depending 
on  their  location  in  phase  space2. 


‘Work  performed  for  USDoE  under  contract  FG03-86ER220 
W.S.  Lawson,  “PDWl  User’s  Manual”,  ERL  Report  April  1984. 
"A.  Fnedman,  et  al ,  APS-DPP,  November  1987. 
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Direct  Implicit  Particle  Simulation  of  Simple  Mirrors.* 
R.  J.  PR0CASSIN1,  J.  C.  CUMMINGS  and  C.  K.  BIRDSALL,  Umver- 
sity  of  California,  Berkeley,  and  B.  I.  COHEN,  Lawrence  Livermore  Na¬ 
tional  Laboratory  —  The  simulation  of  particle  confinement  in  a  simple 
mirror  plasma  is  possible  with  the  one-dimensional  (axial),  relativistic, 
electrostatic  particle  code  TESS.  Our  recent  advances  in  the  develop¬ 
ment  of  TESS  include  the  addition  of  a  quasilinear  RF  energy  diffusion 
model  for  both  ICRH  and  relativistic  ECRH  heating  ,  and  the  addition 
of  a  neutral  beam  injection  model  which  simulates  plasma  heating  via 
charge  exchange  and  impact  ionization  processes.  The  results  of  a  self- 
consistent  TESS  simulation  of  particle  confinement  and  transport  in  an 
auxiliary  heated  simple  mirror  plasma  will  be  presented.  These  results 
will  be  compared  to  those  from  a  Monte  Carlo  RF  heating  code  and  a 
Fokker-Planck  code  vis  a  vis  velocity  distribution  function  resolution, 
computational  cost  effectiveness  and  ease  of  operation. 

‘Work  performed  for  USDoE  under  contract  W-7405-ENG-48. 


The  Effect  of  Thermionic  and  Secondary  Electron 
Emission  at  the  Collector  on  Potential  Drop  and  Trans¬ 
port  Through  the  Plasma-Sheath  Region.*  L.  A.  SCHWAGER, 
U.  C.  Berkeley  —  The  region  between  a  Maxwellian  plasma  source  and 
an  absorbing  surface  which  emits  thermionic  or  secondary  electrons  is 
modeled  numerically  with  electrostatic  particle  simulation  and  analyti¬ 
cally  with  a  kinetic  plasma-sheath  equation.  Results  from  these  analyt¬ 
ical  and  numerical  models  for  collector  potential  and  plasma  transport 
show  excellent  agreement.  Our  energy  transport  results  at  critical  emis¬ 
sion  (field  reversal)  and  beyond  compare  well  with  such  works  as  Hobbs 
and  Wesson1  who  assume  r  <  1  (where  r  =  Ti/Tt,  the  ion/electron  source 
temperature  ratio).  Our  potential  drop  across  the  collector  sheath  is 
slightly  more  negative  (by  5%  for  r  =  0.1  and  by  12%  for  r  =  1)  than 
theirs.  The  plasma  source,  with  a  zero  field  boundary,  emits  equal  fluxes 
of  half-Maxwellian  ions  and  electrons  with  specified  mass  and  temper¬ 
ature  ratios.  The  floating  collector  emits  thermionic  or  secondary  elec¬ 
trons,  with  temperatures  of  1%  Te,  which  are  accelerated  into  and  inter¬ 
act  with  the  plasma.  Spatial  profiles  of  potential  and  particle  and  energy 
fluxes  are  shown  from  simulation  and  compared  with  theory. 


‘Work  performed  for  USDoE  under  contract  DE-FG08-86ER53220. 
1G.  D.  Hobbs  and  J.  A.  Wesson,  Plasma  Physics  9,  85,  (1967). 
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